function [pX,pY,OD] = ProcessImageGUI_v1(handles)

% [pX,pY,OD] = ProcessImage(hdfFilePath,autoROI,userROI,doFitData)
% Finds Gaussian-shaped signals on fluorescence or absorption images.
% Returns:
% pX - The parameters for a Gaussian fit to the data in the x (vertical)
% direction.
% pX(1) = Gaussian amplitude
% pX(2) = Peak position (in units of image pixels)
% pX(3) = Standard deviation (in units of image pixels)
% pX(4) = Offset
% pY - Same description as pX
% OD - The image (OD if it's an absorption image, raw image for a fluor.).
cons = Constants();

y = double(h5read(handles.dpath,'/Image/image'));

%%% Should read all this from the HDF file: %%%
autoROI = handles.autoROI;
doFitData = 1;
userROI = handles.ROI;
handles.img_time = 100e-6;
imgConfigStruct = ParamImgSys(GetImgConfig(1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dim = size(y);
abs = 0;

% Filter image to find ROI
N = 60;
sig = 30;

if size(dim,2) == 3 % Absorption
    num = y(:,:,1) - y(:,:,3);
    den = y(:,:,2) - y(:,:,3);
    den(den==0) = 1;
    ODlin = (num./den)';
    
    OD = -log(ODlin);
    OD(isinf(OD)) = 0;
    OD(isnan(OD)) = 0;
    OD = real(OD);
    OD = OD - mean(mean(OD));
    isAbsImg = true;
elseif size(dim,2) == 2 % Fluorescence
    OD = y';
    OD = OD - mean(mean(OD));
    isAbsImg = false;
else
    error('Can''t determine imaging method');
end

ODfilt = conv2(OD,gaussian2d(N,sig),'same');

yC = 1:imgConfigStruct.short_px_nb;
xC = 1:imgConfigStruct.long_px_nb;

autoROIok = false;

if autoROI
    % Find ROI and extract first approx. signal from filtered image.
    mx = max(max(ODfilt));
    [I,J] = find(ODfilt == mx);
    wndw = 50;
    
    if ~(I - wndw < 1 || I + wndw > max(xC) || J - wndw < 1 || J + wndw > max(yC))
        autoROIok = true;
        
        ROI(1) = I(1);
        ROI(2) = J(1);
        
        yIntFilt = sum(ODfilt(:,ROI(2)-wndw:ROI(2)+wndw),2);
        xIntFilt = sum(ODfilt(ROI(1)-wndw:ROI(1)+wndw,:),1);
        
        ampGuessY = max(yIntFilt) - min(yIntFilt);
        ampGuessX = max(xIntFilt) - min(xIntFilt);
        
        options = optimset('Display','off','TolFun',1e-16,'TolX',1e-16);
        pYtmp = fminsearch(@(p)FitGauss(p,yC,yIntFilt),[ampGuessY ROI(1) 100 min(yIntFilt)],options);
        pXtmp = fminsearch(@(p)FitGauss(p,xC,xIntFilt'),[ampGuessX ROI(2) 100 min(xIntFilt)],options);
        
        % Use data from filtered image to extract data from raw image.
        sigmas = 3;
        xWndw = round(sigmas*pXtmp(3));
        yWndw = round(sigmas*pYtmp(3));
        xROI = ROI(2)-xWndw:ROI(2)+xWndw;
        yROI = ROI(1)-yWndw:ROI(1)+yWndw;
        
        % Remove parts of the ROI if it goes outside of the photo.
        if min(xROI) < 1
            xROI = xROI(xROI > 0);
        end
        if min(yROI) < 1
            yROI = yROI(yROI > 0);
        end
        
        if max(yROI) > imgConfigStruct.short_px_nb
            yROI = yROI(yROI < imgConfigStruct.short_px_nb + 1);
        end
        if max(xROI) > imgConfigStruct.long_px_nb
            xROI = xROI(xROI < imgConfigStruct.long_px_nb + 1);
        end
        
        yInt = sum(OD(yROI,xROI),2);
        xInt = sum(OD(yROI,xROI),1);
        
    else
        warndlg('Can''t autoROI: Estimated signal close to edge of image')
        autoROIok = false;
    end
end
    
if ~autoROIok
    % Use user-input ROI parameters to guess fitting parameters
    yROI = userROI(3):userROI(4);
    xROI = userROI(1):userROI(2);
    
    yInt = sum(OD(yROI,xROI),2);
    xInt = sum(OD(yROI,xROI),1);
    
    pYtmp(1) = max(yInt) - min(yInt);
    pXtmp(1) = max(xInt) - min(xInt);
    
    pYtmp(2:3) = [mean(userROI(3:4)) 150];
    pXtmp(2:3) = [mean(userROI(1:2)) 150];
end

if doFitData
    options = optimset('Display','off','TolFun',1e-16,'TolX',1e-16);
    pY = fminsearch(@(p)FitGauss(p,yROI,yInt),[pYtmp(1:3) min(yInt)],options);
    pX = fminsearch(@(p)FitGauss(p,xROI,xInt'),[pXtmp(1:3) min(xInt)],options);
else
    pY = zeros(1,4);
    pX = pY;
end

axes(handles.axes_processed)
% oneMat = ones(size(OD));
% zeroMat = nan(size(OD));
% zeroMat(yROI,xROI) = 1;
% ODwindow = OD(yROI,xROI);
% ODdisp = OD./max(max(OD));
% ODwindow = ODwindow./max(max(ODwindow));
% ODdisp(yROI,xROI) = ODwindow;
imagesc(OD);
set(gca,'FontSize',10)

line([xROI(1) xROI(end)],[yROI(1) yROI(1)])
line([xROI(1) xROI(1)],[yROI(1) yROI(end)])
line([xROI(end) xROI(end)],[yROI(1) yROI(end)])
line([xROI(1) xROI(end)],[yROI(end) yROI(end)])
axis equal
xlims = xlim;
ylims = ylim;

axes(handles.axes_y_proj)
plot(yInt,yROI,gauss(pY,yC),yC);
set(gca,'YDir','reverse');
ylim(ylims);
set(gca,'FontSize',10)

axes(handles.axes_x_proj)
plot(xROI,xInt,xC,gauss(pX,xC));
xlim(xlims);
set(gca,'FontSize',10)

if isAbsImg
    axes(handles.axes_raw_abs_atoms)
    imagesc(y(:,:,1))
    axis off
    
    axes(handles.axes_raw_abs_no_atoms)
    imagesc(y(:,:,2))
    axis off
    
    axes(handles.axes_raw_background)
    imagesc(y(:,:,3))
    axis off
    
    num_atoms = 0;
    
else
    % Estimate number of atoms for fluorescence imaging
    amp = (pX(1) + pY(1))/2; % Take mean amplitude of two Gaussian fits.
    total_num_electrons = amp*sqrt(2*pi)*mean([pX(3) pY(3)]);
    total_photons = total_num_electrons/imgConfigStruct.coll_eff;
    num_atoms = 2*total_photons/(handles.img_time*cons.Gamma); % num_atoms = (num_photons/time) / (scattering rate)
end

set(handles.text_num_atoms,'string',['# Atoms: ' num2str(num_atoms,'%.2e')]);
set(handles.text_x_amplitude,'string',['X amp.: ' num2str(pX(1),'%.2e')]);
set(handles.text_y_amplitude,'string',['Y amp.: ' num2str(pY(1),'%.2e')]);
set(handles.text_x_width,'string',['X width: ' num2str(pX(3),'%.2e')]);
set(handles.text_y_width,'string',['Y width: ' num2str(pY(3),'%.2e')]);
